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INTRODUCTION 


Resent  testing  was  conducted  at  the  U.S.  Army  Armament  Research,  Development  and 
Engineering  Center  (ARDEC),  Picatinny  Arsenal,  New  Jersey  on  a  soft  recovery  system  for  a 
155-mm  projectile  in  an  assembly  as  shown  in  figure  1  and  with  numerical  details  given  in  figure 
2.  Sense  and  destroy  armor,  such  as  the  155-mm  projectile,  have  very  sensitive  components 
and  packaging  intrinsic  to  their  design.  Projectiles  tested  to  evaluate  performance  and  failure 
analysis  using  a  soft  recovery  system  is  very  expensive.  As  of  date,  recovery  was  employed  to 
retrieve  the  projectile  in  such  a  way  that  the  projectile  does  not  exceed  damage  thresholds.  The 
complexity  of  the  projectile  makes  it  vulnerable  to  damage  during  testing  and  may  become  too 
damaged  for  useful  analysis.  The  use  of  soft  recovery  will  enhance  and  enable  the  verification 
of  launch  functionality/performance  of  the  projectile’s  components,  by  measuring  these  variables 
during  the  actual  testing.  This  report  details  a  non-numerical  closed  form  model  and  its 
predictions  for  the  flow  dynamics  in  the  proposed  soft  recovery  work.  Not  only  is  there  a  need 
for  a  soft  recovery  system,  but  also  for  a  computer  program  that  could  accurately  predict  the 
projectile’s  behavior  as  it  is  fired  through  a  gun  assembly,  which  could  reduce  testing  costs  and 
speed  research  and  development  in  ballistics  design  and  implementation. 
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Scat  gun  assembly  (ref.  1) 
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THEORY 


The  formulas  used  in  the  FORTRAN  program  (app  A)  that  simulates  the  scat  recovery 
system  were  taken  from  a  standard  text  by  Anderson  (ref.  2).  A  shockwave  is  initially  created  as 
the  projectile  enters  the  gun  tube  with  a  velocity  of  516  m/s.  The  shockwave  travels  back  and 
forth  along  the  tube  where  the  length  of  the  path  in  a  given  direction  changes  because  of  the 
projectile’s  forward  motion.  The  Mach  number  of  the  shockwave  changes  from  its  initial  value  to 
when  the  shockwave  reflects  off  the  back  of  the  diaphragm  and  the  front  of  the  piston.  When  the 
pressure  due  to  the  shockwave  exceeds  the  threshold  of  the  diaphragm,  the  diaphragm  bursts 
and  pressure  is  then  exerted  upon  the  back  of  the  piston.  Shockwaves  are  reflected  off  the  back 
of  the  piston  and  the  front  of  the  projectile  until  the  piston  displacement  creates  pressure  on  the 
water  plug  until  it  bursts.  The  energy  expended  by  the  piston  in  bursting  the  water  plug  is  equal 
work  done  in  bursting  the  water  plugs  encasing  the  water  section  of  the  scat  assembly  and 
pushing  the  water  mass  forward.  The  yield  strength  of  the  water  plug  is  used  to  calculate  the 
energy  expended  and  lost  by  the  piston.  The  yield  strength  of  the  diaphragm  and  water  plugs 
were  both  assumed  to  be  equal  to  3,447,378  P.  The  FORTRAN  program  used  in  the 
simulations  has  two  loops:  a  G  loop  with  and  embedded  N  loop.  Both  loops  have  an  initial  value 
of  one  and  are  incremented  by  one  for  each  loop.  Different  Gs  represent  changing  conditions 
when  a  barrier  is  burst.  Each  N  calculates  changing  variables  (i.e.,  Mach  number,  shockwave 
velocity,  projectile/piston  displacements  and  velocities)  during  the  simulation  when  the 
shockwave  completes  a  circuit  in  both  forward  and  backward  directions  as  it  reflects  off  of 
surfaces  within  the  assembly. 

Calculation  of  Initial  Mach  Number 

The  non-reflected  Mach  number  can  be  calculated  by  knowing  the  particle  velocity  in 
front  of  the  projectile.  The  particle  velocity  can  be  approximated  as  the  initial  speed  of  the 
projectile  Vp.  The  initial  Mach  number  is  calculated  as  (fig.  3) 


=  ^siP2  Pi -1)^1 /{r + 0^{P2  ^Pi+ir-  ■/)/(r + ^))  (1) 

where  y  =  c^,  /c„,  K^is  the  shockwave  velocity,  p2  \s  the  shockwave  pressure  in  back  of  the 
shockwave  and  on  the  front  of  the  projectile,  and  p/  is  the  pressure  in  front  of  the  shockwave 
and  on  the  back  of  the  diaphragm. 

The  initial  shockwave  Mach  number  is  given  by 

M -  yl{{/  +  /)/ 2/){p2  /  Pi  -1)  +  1  (2) 

Solving  equation  1  for  p2  /  pi 


{P2/ Pi? -P2/ Pi[2+{yp/Vsf\7  +  i)/r)-{r-i)/r{yp/Vsf  =0 
{P2  /  PI  ?  -P2/Pi(2  +  {Vp  /  Vs  ?{r  +  i)/r)-{r  -  i)/r{Vp  /Vs?=  o 
c>  =  i,b  =  -(2  +  {Vp/Vs?{r  +  i)/r]c  =  -{r-i)/r{Vp/Vs? 
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p1<p2  Shockwave  travels  to  back  of  diaphragm 


Figure  3 

Shockwave  propagation  and  pressure/temperature  changes 
Mathematically,  the  solution  to  this  quadratic  equation  is 

{p2  / pi)±  =  {l / 2a'^-h±'J -4ac^  (3) 

Physically,  the  only  possible  solutions  occur  if 

1.  h<0,'Jh^  -4ac  <-b,p2/pi  =(7/2o)^-/>- V/)'^  -4ac^ 

2.  h>0,-4ac>0,‘^h^ -4ac  >  ahs{-h),  P2  /  Pi  ={l  /  h  +  4h^ -4ac^ 

Substitute  equation  3  into  equation  2  to  find  initial  Mach  number. 

Calculation  of  Sound  Speed  in  Region  Ahead  of  Shockwave 

Temperatures  in  regions  in  front  of  and  in  back  of  the  shockwave  are  function  of 
pressures  in  those  respective  regions,  which  is  given  generally  as 

T^2N  =  {'^NPN  /  +  ((/  +  l)^{y-J))P2N  /  +  +  ^)+  A’JA/  ^  Pn)  H) 

when  shockwaves  are  reflected  in  the  forward  direction  and 

'^2N+I  ={'^2NP2NP2N+1^^  +  ((T'  +  f  )/(/ “  ^)V2A^+/  )/((^  +  “  ^)+  P2.V+/  / /^2/V  )  (5) 
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when  shockwaves  are  reflected  in  the  backward  direction.  Temperatures  with  the  higher 
subscripts  are  in  regions  behind  the  shockwave. 

Substitute  the  previous  relation  for  p^/p/^rom  equation  3  to  find  The  initial 

temperature  7)  in  the  region  where  the  shockwave  is  propagating  is  300K.  The  initial  sound 
speed  in  front  of  the  shockwave  can  then  be  calculated  from  Vsi=4Wr 

Calculation  of  the  Time  for  Shockwave  to  T ransverse  Length  of  Tube  in  the  Forward  Direction 

The  time  zlryv  upon  hitting  the  back  of  a  barrier  the  nth  time  is  zlryv  =7-/^5,  where  Lis 
the  initial  distance  from  the  projectile  to  the  barrier  and  is  the  shockwave  velocity. 

Caicuiation  of  How  Far  the  Projectiie  Traveis  within  Tube  When  Shockwave  Reflects  Off  Front 
of  Projectile 

The  projectile  velocity  and  displacement  at  any  time  can  be  taken  from  standard 
formulas  in  physics.  During  the  time  a  shockwave  is  traveling  through  a  medium,  the  pressure 
on  the  front  of  the  projectile  is  constant  and  is  decelerating  it.  The  nth  distance  traveled  during 
time  zlr„  in  the  forward  direction 


”  2 
Ax^=  I  (y{xt^)-p2f^At/ Mp)dt  =  v{xi^){ATf^)- p2NAATt^  / 2Mp  = 

0 

(6) 

-P2NA{L/Vsf  /2Mp 


where  A  is  the  cross-section  area  of  the  projectile,  v(A:jY)is  the  initial  velocity  of  the  projectile 
when  shockwave  is  either  reflected  or  created  at  the  surface  of  the  projectile,  L  is  the  distance 
the  shockwave  has  to  travel  before  hitting  a  barrier,  is  the  shockwave  velocity,  and  Mp\s  the 
projectile  mass. 

Calculation  of  Time  for  Reflected  Shockwave  to  Hit  Front  of  Projectiie 

As  the  shockwave  reflects  off  the  diaphragm,  the  projectile  is  still  moving  forward,  but  in 
a  decelerating  manner.  The  time  the  shockwave  takes  to  reflect  off  the  diaphragm  and  hit  the 
front  part  of  the  projectile  can  easily  be  shown  to  be 

+  I  ~ ^  BNpVs  (^) 

and  how  far  the  projectile  is  displaced  during  this  time  interval  is 

+  Z  P2N'^'  ^  =  P2N'^{^'’^N  +  l)  /2Mp 

0 


-K'^A'+/)(7'-'4vyv  -zlYyv^/)/ M  -  p  2  A{{L  -  Ax  -Axj^^i)/ M  pfjpVpY  /  2Mp 
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where  v(A:yv+/)is  the  projectile’s  velocity  when  the  shockwave  is  reflected  backwards  and 
is  the  Mach  number  calculated  from  equation  9.  Solving  Axf^^i  involves  a  quadratic  equation. 

^NPIN^P  +  H^/V+/)/ ^BNR^S  ~  P2N^i^~  p(^BNR^s)^)~ 

^BNR^S'^  P2N^P^  2  M  p{{L- Axj^)/ =0 

3  =  P2N^^{^^pi^BNR^s)^) 

b=  {l  +  y{xnj+i)/ Mp}^pVs-p2f^A{L-Axf^)/{Mp{Mgf^pysY^ 

C  =  -v{xfj^.i\L- Axf])/ Mgf^pVs  +  P2N^\^^p{{^~^n)^  ^BNR^s)^) 

^xi\l+]{±)  =  ^-b±'^h^  -4ac'j/{2a)  (8) 

Physically,  the  only  possible  solutions  occur  if 

1.  b<0,^lh^  -4ac  <  -h,  P2  /  Pi  =  V / 2a^-h-ylh^  -4ac^ 

2.  b<0,  —  4ac  >  0,  "Jb^  -4ac  >  abs{-b),  P2  Pi  —  ^  2a^-b  +  'Jh^  -4ac^ 

Knowing  Ar^yand  Axpj_j.  zlr^y+ycan  be  determined  from  equation  8. 

Calculation  of  the  Reflective  Mach  Number  yW/f  after  Shockwave  Hits  Surface  of  a  Barrier 

Upon  forward  reflection  from  a  surface,  the  pressure  in  back  of  the  shockwave  becomes 
P2N+1 3nd  in  the  front  of  the  shockwave  p^^s  shown  in  figure  2.  For  backward  reflection  from 
a  surface,  the  pressure  in  back  of  the  shockwave  becomes  p^/v+z  3nd  in  the  front  of  the 
shockwave  pf^ . 

One  can  now  calculate  the  new  reflected  Mach  number  using 

M^Mc-Mp-Mc^O  (9) 

where  Me  =  (M/  /[m^  -  +  2{y  -  l){Mf  +  A// )/(;'  + 1)^  and  A//  is  the  old  Mach  number 

of  the  shockwave  before  reflection  and  Mp  \s  the  new  reflected  Mach  number. 
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Generally  speaking,  knowing  ,  the  new  pressure  on  the  surface  of  a  barrier  can  be 
estimated.  For  backward  reflection,  F2/V+/  -  FACTOR* P2i^ :\o\  reflection 

-  FACTOR  *  P2N=i ,  where  factor  is  a  constant  obtainable  from  standard  tables  (ref.  3).  In 
both  cases,  the  pressure  with  the  higher  subscript  is  the  shockwave  pressure  on  the  surface  of  a 
barrier  nearest  to  the  reflected  shockwave. 

Calculation  Velocity  of  Piston  as  a  Function  of  its  Displacement 

Computer  simulation  showed  that  the  pressure  Pof  the  shockwave  on  the  piston  is  always 
much  greater  than  the  gas  pressure  on  the  front  of  the  piston  and  the  friction  forces  F  during  the 
time  the  shockwave  travels  towards  or  away  from  the  piston.  Thus,  the  acceleration  of  the 
piston  can  be  approximated  as  constant. 

Aps  =  [{PA-F]/  M  PS  (10) 

There  are  three  forces  acting  on  the  piston:  shock  pressure,  friction,  and  pressure  from 
the  ideal  gas  in  the  chamber  between  the  piston  and  the  water  plug.  The  net  force  on  the  piston 
can  be  used  to  calculate  the  increase  in  its  kinetic  energy  when  the  net  force  on  the  piston 
increases  its  kinetic  energy  as  it  is  displaced  a  distance  X . 

X  XX 

MpsV^S^2=  l((P-vRT)/{Lpyy-X))A-F)dx=  \{-vRT  /{Lpfy  -  X))Adx+  j(PA-  F)dx 
0  0  0 

PMpsV^S  /  ^  =  vRTALN(]  -X/  Lpw)  +  {PA  -  F)X  (11) 

Vps  =  yJ{2vRTApLN{l -X  / Lpw)  +  2{PAp  - F)X  / Mps 
For  the  constants  used,  X / LPW  -.0736  so  the  above  simplified  to 

Vps=4^{PA^^F)x71d^  ox  X  =  MpsV^s/iAPA-F))  (12) 

After  diaphragm  bursts,  the  piston  is  being  accelerated  by  pressure  of  the  shockwave 
from  the  projectile  and  decelerated  by  friction  forces  and  gas  pressure  between  the  piston  and 
the  water  plug.  When  the  water  seal  bursts,  energy  will  be  expanded  to  break  the  seal  which 
decelerated  the  piston  to  zero.  To  calculate  when  the  water  seal  burst,  one  calculates  how  far 
the  piston  travels  to  the  point  when  the  gas  pressure  in  the  region  in  between  the  piston  and  the 
seal  reached  its  yield  strength  Yp, .  If  the  water  plug  bursts  before  the  piston  is  displaced  .v.  the 
critical  displacement  can  be  calculated  using  the  Ideal  Gas  Law  PV  =  vRT. 

Yp=vRT/{A{Lpw-x,,i,)) 


{ALpw  ~  \ni ))  =  vRT  /  Yp) 


(13) 


.V,.,./,  =  Lpyy  -vRT  /(Yp,A) 
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where  v  equals  the  number  of  moles  of  air  between  the  piston  and  water  plug,  T  =  temperature, 
A  =  the  cross-section  area  of  the  front  of  piston,  Xcru  is  the  critical  displacement  of  the  piston 
when  the  water  plug  bursts,  and  Lpfy  =  initial  distance  between  piston  and  water  plug. 

Time/Displacement  Offsetting  of  the  Projectile  and  Piston 

Time  increments  in  the  computer  program  are  not  calculated  continuously;  they 
are,  in  face  only,  calculated  when  the  shockwave  strikes  a  surface  within  the  assembly;  i.e.,  the 
surface  of  the  projectile,  diaphragm,  piston,  or  water  plug.  When  a  barrier  bursts,  the  computer 
program  calculates,  using  a  subroutine,  when  and  where  the  projectile/piston  is  when  the 
bursting  occurs. 


RESULTS  OF  COMPUTER  SIMULATIONS 

The  computer  simulations  and  data  plots  are  shown  in  figure  4.  There  were  few 
quantitative  agreements  with  the  results  of  the  FORTRAN  computer  simulations  and  the 
experimental  data.  The  results  are,  however,  qualitatively  similar  to  the  point  that  there  was  a 
sharp  spike  in  the  pressure  followed  by  a  plateau  outing. 

The  initial  pressure  on  the  projectile  was  calculated  from  the  FORTRAN  program  (app  A) 
to  be  569,000  P,  which  was  consistent  with  the  value  at  point  2  on  the  data  plot  in  figure  4.  The 
initial  pressure  was  determined  to  be  instantaneous  from  the  computer  simulation;  from  the  data 
plot  from  points  1  and  2,  the  plot  shows  a  sharp  pressure  gradient  0.04  sec  with  noisy  entries, 
which  was  the  result  of  the  projectile  not  entering  the  scat  assembly  before  that  time.  In  the 
computer  simulation,  the  projectile  pressure  remains  constant  until  the  diaphragm  bursts  during 
which  the  shockwave  pressure  “spikes”  upwards  when  the  shockwave  hits  the  back  of  the 
diaphragm.  However,  the  data  projectile  pressure  abruptly  decreases  after  point  2  in  figure  4  to 
one  atmosphere,  which  was  the  projectile  pressure  before  it  entered  the  assembly.  This  could 
only  have  happened  if  the  stiff  coupler  in  figure  1  leaked  out  the  pressure  in  the  chamber  in 
between  the  projectile  and  the  stiff  coupler.  The  diaphragm  bursts  at  0.2  sec  in  the  data  plot  and 
0.38  sec  in  the  computer  simulations.  The  difference  was  partially  due  to  the  effect  of  the  stiff 
coupler  in  the  scat  assembly,  which  was  not  noted  of  the  FORTRAN  program.  Differences  of 
pressure  maximum  at  points  3  and  4  in  figure  4  were  most  likely  due  to  the  effects  of  turbulence, 
which  the  computer  program  didn’t  simulate.  For  both  the  data  and  computer  simulation,  the 
projectile  pressure  remained  constant  -  indicating  that  the  water  plugs  burst  during  the  time  after 
the  shock  reflected  off  the  back  of  the  piston  and  hit  the  projectile.  The  negative  pressures  are 
artifacts  from  the  circuitry  in  the  pressure  sensor  since  negative  pressures  are  not  physically 
realistic.  Most  likely,  they  were  offset  incorrectly  when  the  circuitry  processed  the  data  resulting 
in  a  negative  pressure.  Simulations  predicted  that  the  piston  will  decelerate  to  zero  within  the 
passably  which  was  consistent  with  the  experimental  results. 
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The  pressure  in  back  of  the  shockwave  is  always  greater  than  the  pressure  in  front  of  the 
shockwave.  In  figure  4,  the  diaphragm  bursts  when  the  projectile  pressure  is  just  above  the 
yield  strength  of  the  diaphragm  in  the  simulation  plot  and  about  36%  below  it  in  the  case  of  the 
data  plot.  In  both  plots,  the  diaphragm  could  have  only  burst  when  the  shockwave  hit  the  back 
of  the  diaphragm.  The  large  difference  between  the  two  plots  in  projectile  pressure  between  the 
diaphragm  burst  may  have  been  due  the  effect  of  turbulence  in  increasing  pressure  at  surfaces 
during  shockwave  propagation.  Different  values  of  /  could  explain  the  difference,  but  the  same 

value  of  X  =  1-4  was  used. 


156-mm  Projectile 


Figure  4 

Projectile  pressure  versus  time 


CONCLUSIONS 

The  computer  simulations  matched  the  experimental  data  only  in  a  qualitative  sense  to 
the  point  that  there  was  a  sharp  spike  in  the  projectile  pressure  when  the  diaphragm  bursts, 
which  instantaneously  dropped  to  the  initial  projectile  pressure  until  the  water  plug  burst.  The 
disparity  between  the  two  plots  exists  because  the  simulation  failed  to  take  into  account  the 
effects  of  turbulence  on  the  generation  of  pressures  throughout  the  assembly  so  there  was  not 
an  exact  quantitative  concurrence  between  the  simulations  and  data  plots  in  terms  of  projectile 
pressure  and  the  time  coordinate  when  the  diaphragm  bursts.  The  simulations  accurately 
predict  a  zero  exit  velocity  for  the  piston. 
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APPENDIX  A 

FORTRAN  COMPUTER  PROGRAM 
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A  flowchart  for  the  FORTRAN  program  used  in  this  report  is  shown  in  figure  A-1 


Figure  A-1 
Computer  flowchart 
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PROGRAM  SCAT 


IPROGRAMMER:  KENNETH  P.  WALSH,  PHD 
!DATE  08/10/2008 

IPORGRM  CALCULTES  THE  PRESSURE  ON  THE  FRONT  OF  A 
IPROJECTILE  IN  A  SCAT  GUN  ASEMBLY 

INTEGER  Z,N 

REAL  PENETRATION, VO, VS, MASS, PENETRATION_OLD,B 
REAL  P(100),DT(100),DX(100),MC(100),MR(100) 

REAL  GAMMA,R,T(IOO),DTTOTAL,FACTOR,L,A,C,DP_OLD 
REAL  VSOUND,TINIT,VINIT,DELTA_V,NU,RATIO,VP  OLD 
REAL  YI ELD,AREA,V(  1 00),DTTOTAL_OLD,A1  ,BI  ,C  1 
CHARACTER*8  FILENAME 

UNITIALIZE  PHYSICAL  PARAMETERS 

YIELD  =  3447378  lYIELD  STRENGTH  OF  DIAPHRAGM 

MU_ST_STEEL  =  .74 

MU_KN_STEEL  =  .57 

NU=.04I1 

MP  =  5.01 

MASS=45.6  !MASS  OF  PROJECTILE 
ATM=I01325  lATMOSPHERIC  PRESSURE 
P(1)=ATM 

DIA  THICKNESS  =  .002  IDIAPHRAGM  THICKNESS 

LPW=24.4  IDISTANCE  FROM  PISTON  TO  WATER  PLUG 

AREA=  .0189  !CORSS  SECTION  AREA  OF  PISTON  AND  PROJECTILE 

R=289 

T(l)=300 

GAMMA=1.4 

VP_OLD  =  0 
DP_OLD  =  0 

55  FORMAT(E16.7,EI6.6) 

PRINT  *,'ENTER  NUMBER  OF  SECTIONS' 

READ  *,GFINAL 
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!G  CHANES  WHEN  DIAPHRAM  BURST 

!N  CHANGES  AS  SHOCKWAVE  TRAVELS  FROM  ONE  SURFACE  TO 
ANOTHER 

DO  100  G=1,GFINAL,1 
IF(G.EQ,1)THEN 
PRINT  VENTER  FILENAME* 

READ  *,  FILENAME 

OPEN(UNIT=l,FILE=FILENAME,STATUS='NEW') 

PRINT  *, 'ENTER  INTIAL  PROJECTILE  VELOCITY* 

READ  *,VO 
V(l)  =  VO 
TINIT  =  0 

IF(P(1).GT.YIELD)THEN 
GFINAL  =  GFINAL+1 
ENDIF 
ENDIF 

IF(G.EQ.2)THEN 

DELTA_V  =  SQRT(2*YIELD*AREA*DIA_THICKNESS/MP) 
TINIT_OLD  =  TINIT 
ENDIF 

PRINT  venter  length  OF  SECTION* 

READ  *,L 


N=1 

Z=1 

IPENETRATION  =  HOW  FAR  PROJECTLE  IS  DISPLACED 

PENETRATION_OLD  =0 

PENETRATION=0 

DTTOTAL_OLD=DTTOTAL 

!CALCULATE  VELOCITY  OF  SOUN 

VSOUND=SQRT(GAMMA*R*T(l)) 
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DO  WHILE(N.LT.IOO) 

NOLD=N 

IF(N.EQ.1)THEN 

ICALCULATE  INITIAL  MACH  NUMBER 


Al=l 

Bl=-2  -  (GAMMA+l)*V(N)*V(N)*GAMMA/(2*VSOUND*VSOUND) 

Cl=l-(GAMMA-l)*V(N)*V(N)*GAMMA/(2*VSOUND*VSOUND) 

RATIO=(l/(2*Al))*(-Bl-SQRT(Bl*Bl-4*Al*Cl)) 

MR(N)  =  SQRT((GAMMA+l)*(RATIO-l)/(2*GAMMA)  +1) 

T0T01=  MR(N)/(MR(N)*MR(N)-1) 
T0T02=(1+2*(GAMMA-1))*(MR(N)*MR(N)-1) 

T0T03=  (GAMMA+1/(MR(N)*MR(N)))/((GAMMA+1)*(GAMMA+1)) 

IF(-B1.GT.SQRT(B1*B1-4*A1*C1))THEN 

IF(MR(N).GT.1)THEN 

MC(N)=(TOTO  1  )*SQRT(T0T02*T0T03) 

MR(N+1)=  .5*(1+SQRT(1+  4*MC(N))) 

ENDIF 

IF(MR(N).LE.1)THEN 

RATIO=(l/(2*Al))*(-Bl+SQRT(Bl*Bl-4*Al*Cl)) 

MR(N)  =  SQRT((GAMMA+l)*(RATIO-l)/(2*GAMMA)  +1) 

T0T01=  MR(N)/(MR(N)*MR(N)-1) 
T0T02=(1+2*(GAMMA-1))*(MR(N)*MR(N)-1) 

T0T03=  (GAMMA+1/(MR(N)*MR(N)))/((GAMMA+1)*(GAMMA+1)) 
MC(N)=(TOTO  1  )*SQRT(T0T02*T0T03) 

MR(N+1)=  .5*(1+SQRT(1+  4*MC(N))) 

ENDIF 

ENDIF 

IF(-B  1  .LT.SQRT(B  1  *B  1  -4*  A 1  *C  1  ))THEN 

MC(N)=(T0T01)*SQRT(T0T02*T0T03) 

MR(N+1)=  .5*(1+SQRT(1+  4*MC(N))) 
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ENDIF 


VS=MR(N)*VSOUND 

P(2*N)=P(Z)*(GAMMA+1+2*GAMMA*(MR(Z)*MR(Z)-1))/(GAMIV1A+1) 

IF(G.LT.2)THEN 

WRITE(1,55)TINIT,P(2*N) 

ENDIF 

RATIO=(GAMMA+l+2*GAMMA*SQRT(IVlR(N+l)*MR(N+l)-l))/(GAMIVIA+l) 
Rl=((GAMMA+l)/(GAMMA-l)+RATIO)/(l+(GAMMA+l)*RATIO/(GAIVlMA-l)) 
ALPHA=RATIO*Rl 
T(2*N)=  ALPHA*T(2*N-1) 

ENDIF 

IF(N.GT.1)THEN 

ICALCULATE  REFLECTION  MACH  NUMBER 

T0T01=  MR(Z-1)/(MR(Z-1)*MR(Z-1)-1) 
T0T02=(1+2*(GAMMA-1))*(MR(Z-1)*MR(Z-1)-1) 

T0T03=  (GAMMA+1/(MR(Z-1)*MR(Z-1)))/((GAMMA+1)*(GAMMA+1)) 

MC(Z- 1  )=(TOTO  1  )*SQRT(T0T02*T0T03) 

MR(Z)=  .5*(1+SQRT(1+  4*MC(Z-1))) 

PRINT  *,'MR(Z),MC(Z-I),MR(Z+1)' 

PRINT  *,MR(Z-I),MC(Z-1),MR(Z) 

PRINT  *, 'ENTER  FACTOR  TO  CALCULATE  NEWPRESSURE' 

READ  *,  FACTOR 

P(2*N)=P{2*N-l)*FACTOR 

CALL  SOUNDSPEEDl(MR,N,T,VSOUND) 

VS=VSOUND*MR(Z) 

ENDIF 

IF(N.EQ.I)THEN 

DT(Z)=  (L-PENETRATION)A^S  +  TINIT_OLD 
ELSEIF(N.GT.1)THEN 
DT(Z)=  (L-PENETRATION)/VS 
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ENDI 

DXl=V(Z)*(L-PENETRATION)A^S 

DX2=-P(2*N)*AREA*(L-PENETRATION)*(L-PENETRATION)/(2*IVIASS*VS*VS) 

!CALCULATE  DISPLACEMENT  OF  PROJECTILE  AS  SHOCKWAV  E  MOVES  LEFT 
TO  RIGHT 

DX(Z)  =  DXl  +  DX2 

PENETRATION  =  PENETRATION_OLD  +  DX(Z) 

IF(G.EQ.2.AND.N.GE.2)THEN 
PRINT  *,’G.EQ.2.AND.N.GE.2' 

ICALCULATE  POSITION  OF  PISTON  WHEN  WATER  PLUG  BURSTS 
CALL  PSHIFT2(DCRIT,P,N,VP,DP,VP  OLD,DP_OLD,DT,Z,TC,DP 
IF(DP.GE.DCRIT)THEN 
L  =  L+  DCRIT 

PRINT  *, 'WATER  PLUID  HAS  BURST' 

VINIT  =  V(Z) 

TIN1T=  DTTOTAL  OLD  +  TC  +  TINIT  OLD 
WRITE(1,55)TINIT,P(2*N) 

CALL  PISTON2(N,P,VP,LPW,XCRIT) 

GO  TO  100 

ELSEIF(DP.LT.DCRIT)THEN 
PENETRATION  =  PENETRATION  -  DPP 
ENDIF 
ENDIF 

PENETRATION_OLD=PENETRATION 
17  FORMAT(1X,'DX(',I2,')=',E10.4) 

V(Z+1)=V(Z)-  P(2*N)*AREA*DT(Z)/MASS 


18 


51  FORMAT(lX,’V(',I2,')=  ’,E10.4) 


IF(N.GE.2)THEN 

ICALCULATE  REFLECTIVE  MACH  NUMBER 

T0T01=  MR(Z)/(MR(Z)*MR(Z)-1) 
T0T02=(1+2*(GAMMA-1))*(MR(Z)*MR(Z)-1) 

T0T03=  (GAMMA+1/(MR(Z)*MR(Z)))/((GAMMA+1)*(GAMMA+1)) 

MC(Z)=(TOTO  1  )*SQRT(T0T02*T0T03) 

MR(Z+1)=  .5*(1+SQRT(1+  4*MC(Z))) 

PRINT  *,'MR(Z),MC(Z),MR(Z+1)’ 

PRINT  *,MR(Z),MC(Z),MR(Z+1) 

PRINT  *, 'ENTER  FACTOR  TO  CALCULATE  NEWPRESSURE' 

READ  *,  FACTOR 
P(2*N+l)=P(2*N)*FACTOR 

IF(P(2*N+1).GT.YIELD.AND.G.EQ.1)THEN 

ICALCULATE  POSITION  OF  PROJECTLE  WHEN  DIAPHRAGM  BURSTS 

CALL  XTSHIFT1C(TFINAL,DT,DX,L,V,P,N,Z) 

PRINT  *, 'THERE  IS  TOO  NUCH  PRESSURE!' 

PRINT  300,L-PENETRATION 

VINIT=V(Z+1) 

TINIT=DTTOTAL_OLD  +  TFINAL 
WRITE(1,55)TINIT,P(2*N) 

GOTO  100 

ENDIF 

ICALCULATE  SPEED  OF  SOUND  WHEN  SHOCKWAVE  TRAVELS  RIGHT  TO 
LEFT 


CALL  SOUNDSPEED2(MR,N,T,VSOUND) 
VS=VSOUND*MR(Z+l) 
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ENDIF 


IF(N.EQ.1)THEN 

ICALCULATE  PRESSURE  IN  BACK  OF  SHOCKWAVE 

PRINT  *,'MR(Z),MC(Z),MR(Z+1)’ 

PRINT  *,MR(Z),MC(Z),MR(Z+1) 

PRINT  *, 'ENTER  FACTOR  TO  CALCULATE  NEWPRESSURE' 

READ  *,  FACTOR 
P(2*N+l)=P(2*N)*FACTOR 

IF(P(2*N+1).GT.YIELD.AND.G.EQ.1)THEN 

PRINT  *, 'THERE  IS  TOO  NUCH  PRESSURE!' 

300  FORMAT(lX,'ADD  ',E19.4,lX,'TO  THE  NEXT  SECTION') 

PRINT  300,L-PENETRATION 

ICALCULATE  POSITION  OF  PROJECTILE  WHEN  DIAPHRAGM  BURSTS 

CALL  XTSHIFTI  C(TFINAL,DT,DX,L,V,P,N,Z) 

VINIT  =  V(Z+1) 

TINT=DTTOTAL_OLD  +  TFINAL 
WRITE(1,55)TINIT,P(2*N) 

PRINT  *,'DTTOTAL_OLD,DT(Z)' 

PRINT  *,DTTOTAL_OLD,DT(Z) 

GO  TO  100 

ENDIF 

CALL  SOUNDSPEED2(MR,N,T,VSOUND) 

VS=VSOUND*MR(Z+l) 

ENDIF 

PRINT  *,'AREA,MASS.MR(Z+l),VS,VSOUND' 

PRINT  *, AREA, MASS, MR(Z+1), VS, VSOUN 
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ICALCULATEDISPLACEMENT  OF  PROJECTILE  WHEN  SHOCKWAVE  MOVES 
RIGHT  TO  LEFT 

A=P(2*N)*AREA/(2*MASS*VS*VS) 

B1=(MASS*VS*VS) 

B=1+V(Z)A^S  -P(2*N)*AREA*(L-PENETRATION)/Bl 
C=-V(Z)*(L-PENETRATION)A^S  +P(2*N)*AREA*(L-PENETRAT10N 
«&)*(L-PENETRATION)/(2*MASS*VS*VS) 

DX(Z+1)  =  (1/(2*A))*(-B+SQRT(B*B-4*A*C)) 


PRINT  *,'A,B,C,DX(Z+1)’ 

PRINT  *,A,B,C,DX(Z+1) 

PENETRATION  =  PENETRATION_OLD  +  DX(Z+1) 

VINIT=V(Z+1) 

DT(Z+1)=  (L  -  PENETRATION)A^S 
IF(G.EQ.2)THEN 

!CALCUALTE  POSITION  OF  PISTON  WHEN  WATER  PLUG  BURSTS 

CALL  PSHIFT(DCRIT,P,N,VP,DP,VP_OLD,DP_OLD,DT,Z,TC,DPP) 
PRINT  *,'PSHIFT' 

IF(DP.GE.DCRIT)THEN 

L  =  L+  DCRIT 

PRINT  *, 'WATER  PLUID  HAS  BURST' 

VINIT  =  V(Z+1) 

TINIT=  DTTOTAL_OLD  +  TC  +  TINIT_OLD 
PRINT  *,'DTTOTAL_OLD,TINIT_OLD,TC' 

PRINT  *,DTTOTAL_OLD,TINIT_OLD,TC 
WRITE(1,55)TINIT,P(2*N) 

CALL  PISTONl(N,P,VP,LPW,XCRIT) 

GO  TO  100 

ELSEIF(DP.LT.DCRIT)THEN 
PENETRATION  =  PENETRATION  -  DPP 
ENDIF 
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ENDIF 

DTTOTAL=  DT(Z+1)  +  DT(Z)  +  DTTOTAL_OLD 
DTTOTAL_OLD  =  DTTOTAL 
WRITE(l,55)DTTOTAL,P(2*N) 

PRINT  *, 'DTTOTAL, P(2*N)' 

PRINT  *, DTTOTAL, P(2*N) 

V(Z+2)=V(Z+1)-  P(2*N)*AREA*DT(Z+1)/MASS 

2  FORMAT(lX,'TIME  SHOCKWAVE  #',I2,1X,'HITS  PROJECTILE  =',E10.2) 

3  FORMAT(lX,'PRESSURE  OF  SHOCKWAVE  #  ',12,'=  ',E10.2) 

PRINT  *,'Z,DTTOTAL,PENETRATION,N,P(2*N),V(Z+l)' 

PRINT  *,  Z,DTTOTAL,PENETRATION,N,P(2*N),V(Z+l) 

PENETRATION_OLD=PENETRATION 

Z  =  Z  +  2 
N  =  N  +  1 


PRINT  *,'SHOCKWAVE  VELOCITY  =' 

PRINT  *,VS 

ENDDO 

100  CONTINUE 
END 

SUBROUTINE  SOUNDSPEEDl(VIR,N,T,VSOUND) 
INTEGER  N 

REAL  RATIO,GAMMA,ALPHA,MR(100),T(100) 
REAL  VSOUND,Rl 


R=289 

GAMMA  =1.4 

RATIO=(GAMMA  +  1  +  2*GAMMA*SQRT(MR(N+1)*MR(N+1)-1))/(GAMMA+1) 

Rl=((GAMMA+l)/(GAMMA-l)+RATIO)/(l+(GAMMA+l)*RATIO/(GAMMA-l)) 

ALPHA=RATIO*Rl 

T(2*N)=  ALPHA’*^T(2*N-1) 
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VSOUND  =  SQRT(GAMMA*R*T(2*N-1)) 

RETURN 

END 


SUBROUTINE  SOUNDSPEED2(MR,N,T,VSOUND) 
REAL  RATIO,GAMM A,ALPH A,MR(1 00),T(1 00) 


R=289 

GAMMA  =1.4 

RATIO  =  (GAMMA  +  1  +  2*GAMMA*SQRT(MR(N+1)*MR(N+1)-1))/(GAMMA+I) 

ALPHA=RATIO*((GAMMA+l)/(GAMMA- 

l)+RATIO)/(l+(GAMMA+l)*RATIO/(GAMM 

&A-1)) 

T(2*N+1)=  ALPHA*T(2*N) 

VSOUND=  SQRT(GAMMA*R*T(2*N+1)) 

RETURN 

END 

SUBROUTINE  XTSHIFTIA(TFINAL,DT,PENETRATION_OLD,L,V,P,N,Z) 


INTEGER  Z,N 

RE  AL  Z1  ,Z2,Z3,PENETRATION_OLD,TFIN  AL,M  ASS,  ARE  A 
REAL  V(100),P(100),DT(100),F,  MU_KN_STEEL 

MASS=45.26 
AREA=.0189 
MU_KN_STEEL  =  .57 
F=9.8*MASS*MU_KN_STEEL 

PRINT  *, 'SHIFT  IT  lA' 

Z1  =(P(2*N)*  AREA+F)/MASS 
Z3=L  -  PENETRATION_OLD 
Z2=-V(Z) 


TFINAL=(-Z2+SQRT(Z2*Z2-4*Z1*Z3))/(2*Z1) 
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RETURN 


END 

SUBROUTINE  XTSHIFTl B(TFINAL,DT,PENETRATION_OLD,L,V,P,N,Z) 
INTEGER  Z 

REALZl,Z2,Z3,PENETRATION_OLD,TFINAL,MASS,AREA,L 
RE AL  V(  1 00),P(  1 00),DT(  1 00),MU_KN_STEEL 

MASS=45.26 
AREA=.0189 
MU_KN_STEEL  =  .57 
F=9.8*MASS*MU_KN_STEEL 

PRINT  VSHIFT  IT  IB' 

Z1=(P(2*N)*AREA+F)/(2*MASS) 

Z3=L  -  PENETRATION_OLD 
Z2=-V(Z) 


PRINT  *,'Z2*Z2-4*Z1*Z3' 

PRINT  *,  Z2*Z2-4*Z1*Z3 

TFINAL=(-Z2+SQRT(Z2*Z2-4*Z1*Z3))/(2*Z1) 

RETURN 

END 

SUBROUTINE  XTSHIFT2(TFINAL,PENETRATION_OLD,L,V,P,N,Z) 
INTEGER  Z,N 

REAL  Zl,Z2,Z3,PENETRATION_OLD,L, MASS, AREA 
REAL  V(100),P(100),TFINAL,MU_KN_STEEL 

MASS=45.26 
AREA=.0189 
MU_KN_STEEL  =  .57 
F=9.8*MASS*MU_KN_STEEL 

Z1=(P(2*N)*AREA+F)/(2*MASS) 

Z3=L  -  PENETRATION_OLD 
Z2=-V(Z+1) 

PRINT  *,Z1,Z2,Z3 

TFINAL=(-Z2+SQRT(Z2*Z2-4*Z1*Z3))/(2*Z1) 
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RETURN 


END 

SUBROUTINE  XTSHIFT1C(TFINAL,DT,DX,L,V,P,N,Z) 
INTEGER  Z,N 

REAL  Z1,Z2,Z3,TFINAL,MASS,AREA,  MU_KN_STEEL 
RE AL  V(  1 00),P(  1 00),DT(  1 00),DX(  1 00) 

MASS=45.26 
AREA=.0189 
MU_KN_STEEL  =  .57 
F=9.8*MASS*MU_KN_STEEL 
PRINT  VSHIFT  IT  1C’ 

Z1  =(P(2*N)*  ARE  A+F)/MASS 
Z3=  DX(Z) 

Z2=-V(Z) 


TFINAL=(-Z2+SQRT(Z2*Z2-4*Z1*Z3))/(2*Z1) 

RETURN 

END 


SUBROUTINE  PISTON I(N,P, VP, LPW,XCRIT) 

INTEGER  N 

REAL  P(IOO), YIELD,  DELTA_V,VFINAL 

REAL  PISTONMASS,WATERMASS,  L,MU_KN_STEEL,FWORK 
REAL  VP,LPW,XCRIT 

PISTONMASS  =  5.01 
WATERMASS  =  100 
AREA  =  .0189 
L=.002 

YIELD  =  3447378 

MU_KN_STEEL  =  .57 

F=9.8*PISTONMASS*MU_KN_STEEL 

FWORK=  (P(2*N+1)*AREA  -F)*  (LPW-XCRIT  +  24.4) 

DELTA_V  =SQRT(2*FWORK/PISTONMASS) 

VFINAL  =  VP  -  DELTA  V 
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IF(VFINAL.GT.O)THEN 


PRINT  VFINAL  PISTON  VELOCITY  IS' 

PRINT  *,  VFINAL 
ENDIF 

IF(VFINAL.LE.O)THEN 

PRINT  'PISTON  WILL  STOP  AT' 

PRINT  *,  24.4  +  (LPW-XCRIT)  -PISTONMASS*VP*VP/(2*(P(2* 
&N+1)*AREA  -F)) 

PRINT  *, 'METERS  FROM  END  OF  SCAT' 

ENDIF 

RETURN 

END 

SUBROUTINE  PSHIFT(DCRIT,P,N,VP,DP,VP_OLD,DP_OLD,DT,Z,TC,DPP) 
INTEGER  N,Z 

REAL  YIELD,MU_KN_STEEL,NU,MP,DCRIT,P(100) 

REAL  AREA,LPW,T,F,A,TC,Al,A2,A3,DP_OLD,VP_OLD 
REAL  DPP,VP,DP,DT(100) 


MP  =  5.01 

LPW=24.4 

AREA=  .0189 

YIELD  =  3447378 

MU_KN_STEEL  =  .57 

NU=.04I1 

T=300 

R=289 

DCRIT  =  LPW  -  NU*R*T/(YIELD*AREA) 
A=(P(2*N+1)*AREA  -  F)/MP 
VP  =  VP_OLD  +  A*DT(Z+1) 

DP=  DP_OLD  +  VP_OLD*DT(Z)  +  A*DT(Z)*DT(Z)/2 
DPP=  VP_OLD*DT(Z)  +  A*DT(Z)*DT(Z)/2 
A1  =  A/2 
A2  =  VP  OLD 
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A3  =  -DPP 

TC  =  (-A2+SQRT(A2**2-4*A1*A3))/(2*A1) 

PRINT  *,’TC' 

PRINT  *,TC 
DP_OLD  =  DP 
VP_OLD  =  VP 

RETURN 

END 

SUBROUTINE  PSHIFT2(DCRIT,P,N,VP,DP,VP_OLD,DP_OLD,DT,Z,TC,DPP) 
INTEGER  N,Z 

REAL  YIELD, IVIU_KN_STEEL,NU,MP,DCRIT,P(1 00) 

REAL  AREA,LPW,T,F,A,TC,Al,A2,A3,DP_OLD,VP_OLD 
REAL  DPP,VP,DP,DT(100) 

MP  =  5.01 

LPW=24.4 

AREA=  .0189 

YIELD  =  3447378 

MU_KN_STEEL  =  .57 

NU=.0411 

T=300 

R=289 

DCRIT  =  LPW  -  NU*R*T/(YIELD*AREA) 

A=(P(2*(N-I)+1)*AREA  -  F)/MP 
VP  =  VP_OLD  +  A*DT(Z+1) 

DP=  DP_OLD  +  VP_OLD*DT(Z+l)  +  A*DT(Z+l)*DT(Z+l)/2 
DPP=  VP_OLD*DT(Z+l)  +  A*DT(Z+l)*DT(Z+l)/2 
A1  =  A/2 
A2  =  VP_OLD 
A3  =  -DPP 


TC  =  (-A2+SQRT(A2**2-4*A1*A3))/(2*A1) 
DP_OLD  =  DP 
VP_OLD  =  VP 
PRINT  *,'TC' 

PRINT  *,TC 

RETURN 

END 
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SUBROUTINE  PISTON2(N,P,V  P,LPW,XCRIT) 

INTEGER  N 

REAL  P(IOO), YIELD,  DELTA_V,VFINAL 

REAL  PISTONMASS,WATERMASS,  L,MU_KN_STEEL,FWORK 
REAL  VP,LPW,XCRIT 

PISTONMASS  =  5.01 
WATERMASS  =  100 
AREA  =  .0189 
L=.002 

YIELD  =  3447378 

MU_KN_STEEL  =  .57 

F=9.8*PISTONMASS*MU_KN_STEEL 

FWORK=(P(2*(N-l)+l)*AREA  -F)*  (LPW-XCRIT  +  24.4) 

DELTA_V  =SQRT(2*FWORK/PISTONMASS) 

VFINAL  =  VP  -  DELTA_V 
PRINT  *, 'VFINAL  =' 

PRINT  *, VFINAL 

IF(VFINAL.GT.O)THEN 

PRINT  *, 'FINAL  PISTON  VELOCITY  IS' 

PRINT  *,  VFINAL 

ENDIF 

IF(VFINAL.LE.O)THEN 

PRINT  *, 'PISTON  WILL  STOP  AT' 

PRINT  *, 'PISTON  WILL  STOP  AT' 

PRINT  *,  24.4  +  (LPW-XCRIT)  -PISTONMASS*VP*VP/(2*(P(2* 
&N+1)*AREA  -F)) 

PRINT  *, 'METERS  FROM  END  OF  SCAT' 

ENDIF 

RETURN 

END 
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GLOSSARY 


At^ 


L 

^  NR 

ys 

Axf^ 


Mjc 

yp 

Mint 

T2 

Tj 


y, 

u 

Mp 

Lpw 

P 

F 

yp 

yps 

X 


Time  shockwave  takes  for  its  nth  transversal  along  the  tube 
Initial  distance  between  projectile  and  diaphragm  or  piston 
Mach  number  of  shockwave  after  the  nth  reflection 
Velocity  of  sound  in  tube 

Displacement  of  projectile  during  the  nth  transversal  of  the  shockwave 
along  the  tube 

Displacement  of  piston  where  water  plug  bursts 

Constant  used  to  determine  the  Mach  number  after  the  nth  reflection  of 
shockwave 

Particle  velocity 

Mach  number  before  shockwave  is  reflected  off  of  a  surface 

Temperature  in  back  of  shockwave  before  shockwave  first  hits  back  of 
diaphragm 

Temperature  in  front  of  shockwave  before  shockwave  first  hits  back  of 
diaphragm 

Sound  speed  at  temperature  Tj 
Particle  velocity 
Mass  of  projectile 

Initial  distance  between  piston  and  water  plug 

Pressure  exerted  on  back  of  piston  by  shockwave 

Friction  force 

Velocity  of  projectile 

Velocity  of  piston 

Displacement 
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yo 


'^PS 
M  PS 

A 

t 

^BNR 

Vsi 

T, 


Yield  pressure  of  diaphragm  and  water  plugs 

Projectile  velocity  after  nth  reflection/creation  of  shockwave  on  the 
surface  of  the  projectile 

Projectile  velocity  after  nth  reflection  of  shockwave  on  the  surface  of  the 
piston/diaphragm 

Piston  acceleration 
Piston  mass 

Cross-section  area  of  projectile  and  piston 
Time 

nth  backward  reflected  Mach  number 
Speed  of  sound  in  region  with  temperature 
Projectile  pressure  in  front  of  the  shockwave 
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